# Purpose: Loads monetized mortality damages. Note that this function loads
# outputs from the next step in the analysis `3_valuation`, but for simplicity
# of the code base and since it's mostly needed for in-text summary stats and
# communications purposes, we keep the code here. Like `impacts.R`, it
# relies upon a set folder structure in `data/3_valuation`.

D_INPUT_DEFAULT = glue('{DB}/3_valuation/impact_region')

ADJ_INPUT_DEFAULT = glue('{DB}/3_valuation/inputs')


#' Pulls mortality damages from valued monte carlo projection output.
#'
#' For this to work, damages must exist as a CSV generated by 
#' calculate_damages.py.
#'
#' @param ssp SSP scenario (SSP2-4)
#' @param rcp RCP scenario (rcp45, rcp85)
#' @param iam Economic modelling scenario (low, high)
#' @param input_dir Directory containing extracted impacts (note default)
#' @param regions Regions, can be IRs or aggregated regions. Also accepts:
#' - all: all ~25k impact regions; 
#' - iso: country-level output; 
#' - global: global outputs
#' @param qtile accepts mean and/or quantiles, e.g., c('mean', 'q05', 'q95').
#' @param valuation Valuation assumptions, e.g., "vsl_epa_scaled" Should match 
#' the column in the `calculate_damages.py` output CSV.
#' @param scn output value. (deaths, costs, deathcosts)
#' @param year_list List of years to extract between 2020 and 2099.
#' @param units dollars (2019 USD) or share_gdp
#' @param scale_variable Can be used to scale output, e.g. millions of dollars
#' @param as.DT Outputs data.table rather than dataframe.
#' @param single_col logical. keeps only requested variables from data
#' @param deryugina_scalar logical. loads valuation input that rescales the age group 3 (oldes) life exepctancy according to results from Deryugina et al (2019).
#' @return Dataframe containing monetized output.
get_mortality_damages = function(
    ssp='SSP3',
    rcp='rcp85',
    iam='low',
    scn='deathcosts',
    valuation='vsl_epa_scaled',
    units='dollars',
    qtile='mean',
    regions='all',
    input_dir=D_INPUT_DEFAULT,
    year_list=seq(2020,2099),
    scale_variable=1,
    ren_var=NULL,
    as.DT=FALSE,
    single_col=TRUE,
    deryugina_scalar=FALSE) {

    # Parse inputs to determine list of regions
    region_list = return_region_list(regions)
    resolution_list = check_resolution(region_list)
    scale_func = function(x, scl) (x * scl)

    # Handle deaths + costs.
    if ("deathcosts" %in% c(scn)) {
        scnlist=c('deaths', 'costs')
        add_costs = TRUE
    } else {
        scnlist=c(scn)
        add_costs=FALSE
    }

    # Generate list of variables.
    varlist = c()
    for (s in c(scnlist))
        for (q in c(qtile))
            varlist = c(varlist, glue("monetized_{s}_{valuation}_{q}"))

    suffix <- ifelse(deryugina_scalar, "_deryugina_scalar", "")

    # Open base file.
    df = memo.csv(glue(
        '{input_dir}/damages_IR_{valuation}_{rcp}_{iam}_{ssp}{suffix}.csv'),
            key=c('region', 'year'))
    
    dflist = list()
    all_vars = c('region', 'year', varlist)

    # Extract Impact Regions
    if (!is.null(resolution_list[['ir_level']]))
        dflist[['ir_level']] = df[
            year %in% year_list  &
            region %in% resolution_list[['ir_level']],
            ..all_vars]

    # Collapse aggregated regions.
    if (!is.null(resolution_list[['aggregated']])) {
        aggregates = get_children(resolution_list[['aggregated']])
        for (reg in names(aggregates)) {
            dflist[[reg]] = df[
                # Subset regions/years.
                year %in% year_list &
                region %in% aggregates[[reg]],
                # Extract columns.
                ..all_vars][
                    # Sum.
                    , lapply(.SD, sum), by=year, .SDcols=varlist][
                        # Replace Region.
                        , region := reg]
        }
    }

    df = rbindlist(dflist, use.names=TRUE)

    # Deaths + Costs.
    if (add_costs) {
        for (q in c(qtile)){
            dcv = glue("monetized_deathcosts_{valuation}_{q}")
            dv = glue("monetized_deaths_{valuation}_{q}")
            cv = glue("monetized_costs_{valuation}_{q}")
            df = df[, (dcv) := get(dv) + get(cv)]
            varlist = c(varlist, dcv)
        }
        all_vars = c('region', 'year', varlist)
    }

    # Calculate share GDP if needed.
    if (units=='share_gdp') {
        gdp_df = get_econvar(units='gdp', regions=regions, 
            year_list=year_list, ssp=ssp, iam=iam, as.DT=TRUE)
        df = merge(df, gdp_df, by=c('region', 'year'))
        sharelist = c()
        for (var in varlist) {
            share = str_replace(var, 'monetized', 'sharegdp')
            df = df[, (share) := get(var) / gdp][, -..var]
            sharelist = c(sharelist, share)
        }
        varlist = sharelist
        all_vars = c('region', 'year', varlist)
    }

    # Scale/reorder.
    df = df = df[, (varlist) := lapply(.SD, scale_func, scl=scale_variable),
        .SDcols=varlist][
            , ..all_vars][
                order(region, year)]

    check = (length(c(scn))==1 & length(c(qtile))==1)
    if (single_col & check)
        df = filter_single_column(df, varlist, scn, qtile, valuation)

    # Convert to dataframe or leave as data.table.
    if (!as.DT) df = data.frame(df)

    return(df)
}

filter_single_column = function(df, varlist, scn, qtile, valuation){
    tag = varlist[grepl(glue("{scn}_{valuation}_{qtile}"), varlist)]
    keep = c('region', 'year', tag)
    return(df[, ..keep])
}
